Collective excitations of trapped Fermi or Bose gases 
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A new method is developed to calculate all excitations of trapped gases using hydrodynamics at 
zero temperature for any equation of state fi = fJ,(n) and for any trapping potential. It is shown 
that a natural scalar product can be defined for the mode functions, by which the wave operator 
is hermitian and the mode functions are orthogonal. It is also shown that the Kohn-modes are 
exact for harmonic trapping in hydrodynamic theory. Excitations for fermions are calculated in the 
BCS-BEC transition region using the equation of state of the mean-field Leggett-model for isotrop 
harmonic trap potential. 
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Several experiments on trapped ultracold gases probed 
in the past decade the collective excitations of atomic 
gases. Earlier measurements on bosons Q, Q and more 
recent measurements on fermions 0, near Feshbach 
resonances can be explained rather satisfactorily using 
hydrodynamics at zero temperature. In his seminal pa- 
per 0] Stringari applied first hydrodynamics for trapped 
bosons undergoing Bose-Einstein condensation. His pre- 
dictions were confirmed by experiments Q, 0- Later, 
using the same approach he predicted also the qualita- 
tive behavior of low lying modes for fermions in the whole 
crossover region from a BCS type superfiuid Fermi gas to 
a molecular Bose-Einstein condensation (BEC) [?], |g, • 
Now, several recent theoretical papers appeared in the 
literature 0EJHHE1OII using hydrodynamic the- 
ory to better explain the measurements on the BCS-BEC 
transition. In general, no exact solution to the hydrody- 
namic equations are known, except when the equation of 
state has the polytropic form fJ,(n) oc rf 1 [TH ]. 

The hydrodynamic approach leads to a wave equa- 
tion for the density oscillations. In principle, this wave 
equation can be solved for a single oscillating mode, if 
the boundary conditions for the density oscillations are 
known. Bulgac et al. [H| has written the eigenvalue 
equation for a single mode in such a way that the two 
sides were hermitian, but did not address the question 
of the function space to which all the excitations should 
belong. Here we use a different approach. For a general 
equation of state /i = /x(n) it is usually very difficult to 
prescribe appropriate boundary conditions at the surface 
of the gas. There are a few examples where this problem 
is circumvented using some ansatz on the spatial forms of 
the excitations Here we shall introduce a natural 

scalar product, by which the wave operator itself is her- 
mitian and the boundary conditions can be treated as in 
quantum-mechanical problems: the mode functions are 



square integrable functions. The scalar product we shall 
use automatically ensures particle conservation. By this 
way finding excitation frequencies are relatively easy: the 
task is to calculate matrix-elements of the wave opera- 
tor with the natural scalar product, and then calculate 
the eigenvalues of the resulting matrix. We shall demon- 
strate the whole procedure for the mean-field model of 
the BCS-BEC transition for isotropic trap potential and 
compare our results with that of the scaling ansatz ap- 
proach for the lowest I = monopole mode. 

In hydrodynamic theory for trapped gases at zero tem- 
perature density oscillations are given by the continuity 
equation 

' >n ■ V(nu) - 0, 



dt 

and the Euler-equation 



5 + uVu=— (VP)- 
at mn 



(W) 



(1) 



(2) 



where u is the velocity field, n is the density, t is the 
time, P is the pressure, m is the particle mass, and V is 
the external trapping potential. 

Knowing the equation of state y, = \x(n) at zero tem- 
perature in the corresponding homogeneous system the 
equilibrium density in the trapped case can be deter- 
mined from the local chemical potential 



H = /i(r) = M(n (r)) = /i - V(r), 



(3) 



where fio is the overall constant chemical potential. For 
confining potentials the solution of this equation for pos- 
itive no(r) supplies an equilibrium density, which has a 
finite support with a well defined boundary. Typically 
no decreases to zero by approaching the boundary. Us- 
ing the thermodynamic identity 



dfi 1 fdP\ 



dn n \ dn ) 



(4) 
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valid also at T = the gradient of Eq. (|2J) gives the 
important vector identity 



-A (r)^M=W(r), Mr) ^ 

no(r) \dn 



(5) 



n=no(r) 



Eq. (JjjJ plays a central role in the following in simplifying 
the linearized hydrodynamics. In mechanical equilibrium 
the equilibrium pressure Pq must satisfy 



VPo(r) = -n (r)W(r), 



(6) 



otherwise the right hand side of Eq. Q will not vanish 
for u = (this is the local form of the Archimedes-law 
for a general external potential). Close to equilibrium u 
and Sn(r,t) = n(r,t) — n (v) are small, and P can be 
expanded to first order in Sn as 



P{v,t)=P {v) + 



dP 

dn 



Sn(r, t). 



(7) 



An important restriction for Sn is particle conservation 
J d 3 r Sn(r, t) = 0. Linearizing the continuity equation 
l (X| l and the Euler-equation (0 in Sn and u using the 
local form of the Archimedes-law © and the identity (J5J 
the linearized hydrodynamic equations can be written as 



— + V(n u) = 0, 



du 



A 



Sn 



n^m 

Let us introduce a new field by 



' Mr) 
n (r) 



Sn(r, t) 



(8) 



(9) 



(10) 



where "J has the same support as no, A$ and Sn. From 
now on, we allow complex fields \& (which are more con- 
venient for problems, where angular momentum is con- 
served). Eliminating u from Eqs. © and © a wave 
equation 



d 2 V 



■G*tf = 0, 



can be derived for \&, where G* is given by 

™o(r) _ 



Gvji — 



' Ao(r) 
n (r) 



%(r) 
"o(r) ' 



(11) 



(12) 



The main advantage of the field \& is that its wave oper- 
ator Gp is manifestly hermitian 16] with respect to the 
scalar product 



(*l|*a> = / d 3 r^l 

Jn a (r)>0 



(r)*a(r) 



(13) 



The scalar product l|13(l is trivially the correct one for a 
homogeneous system with periodic boundary conditions. 
It was used for a weakly interacting trapped Bose-gas 
[rTf. where u(n) oc n. The same idea of finding a proper 
scalar product for eigenmodes of a trapped, noninteract- 
ing Bose or Fermi gas at finite temperature using the 
hydrodynamic approach was applied in 18]. A single 
eigenmodc 

9i(r, t) = sinK-f + <£ )*;(r), (14) 

fulfills the eigenvalue equation 

w?*i(r) = 6**i(r). (15) 

Solutions of l|15fl are (or for degenerate eigenvalues can 
be) orthonormalized with the scalar product i|13|) 



n (r)>0 



d 3 r*|(r)^-(r). 



(16) 



^oO") = Const • y no(r) /Ao(r) is always a solution to 
lfT5Jl with cj = 0. Eq. ((TQl implies that 6n j; = *o*i) 
thus the orthogonality relation (|16fl shows that all the 
modes with i ^ are automatically particle conserving, 
and the mode should be cancelled from the solutions. 

Taking a complete orthonormal basis, i.e., 5i j = 
(<Pi\<Pj) the squared excitation frequencies lo 2 can be ob- 
tained from the eigenvalues of the matrix 



Gi,j = (iPi\Gyi\<Pj). 



(17) 



The matrix elements in (|17|l require the knowledge or the 
numerical evaluation of spatial derivatives of the basis 
functions. Usually this causes big numerical errors be- 
cause the high lying modes are rapidly oscillating func- 
tions. In practice, it is much better to apply the spa- 
tial derivatives to the (spatially varying) coeffectients of 
the wave equation, which are usually not oscillating too 
much. 

The wave operator (|12|l has the structure G$ = 
— RVQVR. If there exists a similar system for which 
the boundary is the same and the wave operator has also 
the structure G° = — i?oVQoVi?o but with known spec- 
tra and eigenfunctions 



gV 



(o) 



(18) 



then one can eliminate the unwanted spatial derivatives 
of the basis functions in the matrix elements if the basis 
is given by (fi, (i — 0, 1, . . .). Let us introduce a and (3 
by 

a = a(r) = Q/Q , P = (3(r) = R/R , (19) 
then the matrix elements can be written as 



G itj = I d- i rip*{r)ip j (r)G ttj (r), 



(20) 
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where 



M , e (o) 



■af3 2 +i?gQ a(V/3)(V/3) 



--^V [Qo/3 2 (Va) 



(21) 



For an isotrop harmonic trapping potential V(r) = 
mu> 2 r 2 /2 and for any equation of state /i = [i(n) a whole 
series of exact solutions of the wave equation can be 
given. If ^(r) is chosen to be 



*(r) = Const a 



np(r) 



(22) 
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then this mode function fulfills the wave equation with 
eigenvalue oj 2 — tu 2 l. The three I — 1 modes are the 
Kohn- modes (see Ref. ^3|) for isotropic trapping. 

As a specific, nontrivial model let us consider the 
mean-field model of Leggett Q for the BCS-BEC transi- 
tion. The Leggett model is fixed in homogeneous systems 
by the gap equation 



^ 2 \E k e k 



and by the number equation 



Anh 2 a ' 



(23) 



(24) 



where E v = ^ '(e k - f i) 2 + A 2 , e k = h 2 k 2 /{2m), A is 
the pairing gap and a is the s-wave scattering length. 
The equation of state /x = fJ,(n) is implicitly given by 
the model. This model captures the essential features of 
the BCS-BEC transition. However, recent Monte-Carlo 
data on the equation of state show that there are 
corrections to the mean-field results of the Leggett model 
which should be taken into account for the equation of 
state, especially close to unitarity (i.e., around the a = oo 
point). Here we study the above model for simplicity. In 
the trapped case we use the n(n) function taken from 
the model and solve for the density profile keeping 
N = J d 3 ru(r) to be fixed. The equilibrium pressure for 
any trap potential V(r) can be calculated from the local 
form of the Archimedes law ©. Knowing the pressure 
and the density the calculation of A$ (r) with help of (JSJ 
is straightforward. The details of the full calculation for 
the Leggett model will be published elsewhere [lflj . 

For isotrop harmonic trapping there is a dimensionless 
coupling parameter: k = d/ (aN 1 / 6 ), where d = y^h/mujo 
is the oscillator length. The spectra depends only on k. 
In three cases the spectra is exactly known 0, 0, |nj 
because the equation of state has a polytropic form: /i oc 
n 1 . These particular values are n = —oo (BCS limit), 
k = (unitarity limit) and k = oo (BEC limit). In 
these cases all the mode functions can be constructed 



FIG. 1: Excitation frequencies = ui/uio of the Leggett model 
as a function of k = d/(aN 1/(i ) on both sides of the BCS-BEC 
transition for isotrop harmonic trapping. Arrows on the left 
(right) denote the excitations for a non-interacting Fermi gas 
(for a weakly interacting Bose gas). 




FIG. 2: The lowest I = excitation frequency Q = ui/uo as 
a function of K = d/(aN 1/6 ) on both sides of the BCS-BEC 
transition for isotrop harmonic trapping. Arrow on the left 
(right) denotes the excitation for a non-interacting Fermi gas 
(for a weakly interacting Bose gas). 



exactly lid . even in the nonisotropic case (The methods 
of Refs. |l7ll2fj| can be easily employed to the polytropic 
equation of state). We used the k — — oo mode functions 
21 1 on the BCS-side and the n = oo mode functions 
11 1 on the BEC-side as basis functions. Our numerical 
results for different angular momentum I can be seen on 
Fig. ^ Arrows on both sides show the limiting well- 
known collective oscillation frequencies 0. In Fig. [5] the 
behavior of the lowest I = quadrupole mode can be seen 
as a function of k. This mode is the lowest k dependent 
mode on Fig. ^ The scaling ansatz approach 0, 
gives quite a good result for this particular mode. On 
the scale of Fig. [21 the two curves would be practically 
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0.0005 



g -0.0005 




the monopole mode and showed that the difference be- 
tween the two frequencies are extremely small in the 
whole region of the BCS-BEC transition for the mean- 
field Leggett model. 
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FIG. 3: Comparison of the exact frequencies and those of 
given by the scaling ansatz as a function of k = d/(aN 1 ^ 6 ) 
for the lowest 1 = monopol mode. 8Q, 2 defined as SQ 2 = 
(lu 2 — uj 2 c )/u>o (See text for ui 2 c ). 

indistinguishable. 

In order to give a quantitative measure about the qual- 
ity of the latter approach we compare in Fig. [3] our 
excitation frequencies with those given by the scaling 
ansatz. In the isotrop harmonic trapping potential case 
o& is given [HE3 by u4/ w o = 9{ndn/dn) / (2(V)) - 1, 
where the average of a quantity like V is taken as (V) — 
J d 3 r V(r)no(r). From the figure it is clearly seen that 
scaling ansatz is exact at k, = 0, ±oo, but between these 
values it is not. However, the difference in the isotropic 
case is so small that is much less than the experimen- 
tal resolution. We have preliminary data [ijj for the 
excitation frequencies for the experimentally relevant ax- 
ially symmetric harmonic trapping as well. Once again 
the scaling ansatz differs a little for the radial and ax- 
ial quadrupole modes for a general intermediate coupling 
k. For such large anisotropics as in 0, 0] however, the 
deviation <5f2 2 is much bigger than in the isotropic case. 

Finally, let us turn to the conclusions. We gave a 
straightforward method how to solve the hydrodynamic 
equations in the trapped case if the equation of state 
is known. We introduced a natural scalar product for 
the transformed wave operator by which the operator is 
hermitian. Collective excitations by our method can be 
found by simply diagonalizing the matrix of the wave op- 
erator on some basis. The power of our method lies in the 
fact that we calculate the whole spectra. We can predict 
the behavior of the excitations also for those modes for 
which no scaling ansatz is known. The method is not lim- 
ited to a particular trap potential (isotropic or not), nor 
a given mode. There is no additional approximation, the 
method calculates the (numerically) exact modes given 
by the hydrodynamic theory. Here we wanted to present 
the basic formalism and applied to the simplest isotrop 
harmonic case. We compared the numerically exact ex- 
citation frequencies with that of the scaling ansatz for 
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